/*!  \file  FullSVD.cpp
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      04/29/2013 05:33:43 PM
 *   \copyright GNU Public License.
 */

#include <iostream>
#include <assert.h>
#include <stdlib.h>
#include <iomanip>
#include <string.h> // memcpy
#include <math.h>
#include <stdio.h>

#include "Matrix.h"
#include "FullSVD.h"
#include "Timer.h"
#include "MatrixGenerator.h"
#include "UnitTester.h"
#include "TriEigen.h"
#include "BiSVD.h"
#include "FormSVD.h"
#include "BiReduction.h"
#include "commonSVD.h"
#include "ThreadPool.h"
using namespace std;

FullSVD::FullSVD(const char * name_): Implementation(name_,FullSVDGroup),
	reduction(NULL), biSVD(NULL),	form(NULL),	
	sanity(false), aColMajor(false), aRowMajor(false) {
		setCompileTime(__DATE__, __TIME__);
	canBeParallel = false;
};
void FullSVD::saveCompute(Matrix & A, Matrix & U, Matrix & S, Matrix & VT){
	reset();
	assert(aColMajor || aRowMajor);
	success &= acceptRowMajor() || acceptColMajor(); assert(success);
	bool dimentionMatch = checkSVDDimentions(A, U, S, VT);
	bool haveMem = checkSVDMemory(A, U, S, VT);
	success &= haveMem && dimentionMatch; assert(success);
	if (!success) return;
	if ((A.isRowMajor()&&acceptRowMajor())||(A.isColMajor()&&acceptColMajor())){
//		cout << " HHHHAAAA " << endl;
//		cout << A.rows() << "x" << A.cols() << " - ";
//		cout << U.rows() << "x" << U.cols() << " * ";
//		cout << S.rows() << "x" << S.cols() << " * ";
//		cout << VT.rows() << "x" << VT.cols() << endl;

		if (!A.sameOrder(U) )  U.transpose();
		if (!A.sameOrder(S) )  S.transpose();
		if (!A.sameOrder(VT)) VT.transpose();

		compute(A, U, S, VT);
		A.transpose();
//		S.transpose();
//		VT.transpose();
//		cout << A.rows() << "x" << A.cols() << " - ";
//		cout << U.rows() << "x" << U.cols() << " * ";
///	cout << S.rows() << "x" << S.cols() << " * ";
//		cout << VT.rows() << "x" << VT.cols() << endl;
	} else {
		A.transpose(); 
		if (!A.sameOrder(U) )  U.transpose();
		if (!A.sameOrder(S) )  S.transpose();
		if (!A.sameOrder(VT)) VT.transpose();
		compute(A, VT, S, U);
		U.transpose(); VT.transpose(); 
		S.transpose();	A.transpose();
	}
	assert(checkSVDDimentions(A, U, S, VT));
};

FullSVD::~FullSVD(){
	if (reduction) delete reduction; reduction = NULL;
	if (biSVD) delete biSVD; biSVD = NULL;
	if (form) delete form; form = NULL;
};

FullSVD::FullSVD(BiReduction * red_, BiSVD * biSVD_, FormSVD * form_): Implementation("no name",FullSVDGroup){
	sanity    = false;
	aColMajor = true;
	aRowMajor = true;
	reduction = red_  ; assert(reduction != NULL);
	biSVD     = biSVD_; assert(biSVD     != NULL);
	form      = form_ ; assert(form      != NULL); 	
	setCompileTime(__DATE__, __TIME__);
	setCompileTime(reduction->getCompileTime());
	setCompileTime(    biSVD->getCompileTime());
	setCompileTime(     form->getCompileTime());
	int l = strlen(reduction->getName());
	l   +=  strlen(biSVD->getName());
	l   +=  strlen(form->getName());
	l   += 10;
	char buf[l];
	sprintf(buf, "%s + %s + %s", reduction->getName(), biSVD->getName(),
			 form->getName() );
	setName(buf);
	canBeParallel = reduction->getCanBeParallel();
	canBeParallel |=    biSVD->getCanBeParallel();
	canBeParallel |=     form->getCanBeParallel();
}

bool FullSVD::setNumberOfThreads(const int nThreads_) {
	bool res = true;
	if (reduction->getCanBeParallel())
		res &= reduction->setNumberOfThreads(nThreads_);
	if (biSVD->getCanBeParallel())
		res &=     biSVD->setNumberOfThreads(nThreads_);
	if (form->getCanBeParallel())
		res &=      form->setNumberOfThreads(nThreads_);
	return res;
}
void FullSVD::compute(Matrix & A, Matrix & U, Matrix & S, Matrix & V ) {
	assert(reduction != NULL && biSVD != NULL && form != NULL);

	assert(A.isGeneral()    && S.isDiagonal() &&
			 U.isOrthogonal() && V.isOrthogonal());

	Matrix Q   (A.rows(), A.rows(), ElementaryReflectors);
	Matrix P   (A.cols(), A.cols(), ElementaryReflectors);
	Matrix B   (A.rows(), A.cols(), BiDiagonal );

	reduction->compute(A, Q, B, P ); // A is destroyed in this computation
	biSVD->saveCompute(B, U, S, V ); // SVD of BiDiagonal matrix
	form->compute     (Q, P, U, V ); // assembly the final result

}


UNITTEST(NormSVDTest);
NormSVDTest::NormSVDTest() {
	setName 				("Test SVD norm computatuon");
	setDescription  	("Test if svd norms are compiuted correctly");
	setTestAuthor		("Yordan Madzhunkov");
	setDeveloperAuthor("Yordan Madzhunkov");
	setCompileTime    (__DATE__, __TIME__ ); 
};
bool NormSVDTest::test() const { 
	// This function make sure that everything that is need to determine
	// if a svd implementaion is valid is implemented correctly
	const int m = 5; const int n = 6;
	double a[] = {8.79,  6.11, -9.15,  9.57, -3.49,  9.84,
	              9.93,  6.91, -7.93,  1.64,  4.02,  0.15,
	              9.83,  5.04,  4.86,  8.83,  9.80, -8.99,
	              5.45, -0.27,  4.85,  0.74, 10.00, -6.02,
	              3.16,  7.98,  3.01,  5.80,  4.27, -5.31};
	double u[] = {
		2.513827927204962e-01,   8.148366860863393e-01,  -2.606185055842206e-01,
		3.967237771305970e-01,  -2.180277636865458e-01,   3.968455517769293e-01,
		3.586615001880026e-01,   7.007682094072525e-01,  -4.507112412166434e-01,
		1.402099498711207e-01,   6.921510074703637e-01,  -2.488880111592842e-01,
	  -2.208114467204375e-01,   2.513211496937533e-01,   5.891194492399430e-01,
		3.661704447722307e-01,  -3.685935379446174e-01,   3.859384831885422e-01,
		4.342486014366712e-01,  -6.265282503648171e-01,   4.076352386533527e-01,
	  -9.796256926688661e-02,  -4.932501428510238e-01,  -6.226840720358041e-01,
	  -4.395516923423331e-01	};
	double s[] = {
	   2.746873241822184e+01,   2.264318500977470e+01,   8.558388228482583e+00, 
		5.985723201512135e+00,   2.014899658715758e+00};
	double v[] = {
	   5.911423764124368e-01,   2.631678147140567e-01,   3.554301738628274e-01, 
		3.142643627269274e-01,   2.299383153647486e-01,   5.507531798028815e-01,
	   3.975667942024255e-01,   2.437990279263305e-01,  -2.223900006854466e-01, 
	  -7.534661509534580e-01,  -3.635896866974969e-01,   1.820347901350359e-01,
	   3.347896906244723e-02,  -6.002725806935829e-01,  -4.508392689223077e-01, 
		2.334496572447148e-01,  -3.054757327479318e-01,   5.361732698764649e-01,
	   4.297069031370181e-01,   2.361668062811257e-01,  -6.858628638738116e-01, 
		3.318600182003101e-01,   1.649276348845101e-01,  -3.896628703606128e-01,
	   4.697479215666589e-01,  -3.508913988837022e-01,   3.874446030996731e-01, 
		1.587355595821562e-01,  -5.182574373535355e-01,  -4.607722286054854e-01,
	  -2.933587584644038e-01,   5.762621191338905e-01,  -2.085291798087080e-02, 
		3.790776670601606e-01,  -6.525516005923975e-01,   1.091068082007294e-01};
	Matrix A    (m,n, General   , a);
	Matrix U_ref(m,m, Orthogonal, u);
	Matrix S_ref(m,n, Diagonal  , s);
	Matrix V_ref(n,n, Orthogonal, v);
	V_ref.transpose();
	const int  nNorms = 6;
	double norm[nNorms];
	const double limit = 5e-14;
	getNorms(norm, A, U_ref, S_ref, V_ref);
	bool res = true; for(int i = 0; i < nNorms; i++) res &= norm[i] < limit;
	return res;
}
FullSVDUnitTest::FullSVDUnitTest(FullSVD & svd_): svd(svd_) {
	int l = strlen(svd.getName());
	char buf[l + 20]; 
	sprintf(buf, "FullSVD %s", svd.getName());
	setName 				(buf);
	setDescription  	("Test registered bi svd implementations");
	setTestAuthor		("Yordan Madzhunkov");
	setDeveloperAuthor("Yordan Madzhunkov");
	setCompileTime(__DATE__, __TIME__); // each class that inherits this one should also have this line.
	setCompileTime(svd.getCompileTime());
};
bool FullSVDUnitTest::test () const {
	const int m = 5; const int n = 6;
	double a[] = {8.79,  6.11, -9.15,  9.57, -3.49,  9.84,
	              9.93,  6.91, -7.93,  1.64,  4.02,  0.15,
	              9.83,  5.04,  4.86,  8.83,  9.80, -8.99,
	              5.45, -0.27,  4.85,  0.74, 10.00, -6.02,
	              3.16,  7.98,  3.01,  5.80,  4.27, -5.31};
	Matrix A (m, n, General, a), U (m, m, Orthogonal);
	Matrix tmp(A,CopyWithData);
	Matrix S (m, n, Diagonal  ), VT(n, n, Orthogonal);
//	tmp.changeOrder();
	svd.saveCompute(tmp, U, S, VT ); // a is destroyed in this computation
//	cout << endl << svd.getName() << endl;
//	cout.precision(16);
//	cout << scientific;
//	cout << S;
	const int  nNorms = 6;
	double norm[nNorms];
	getNorms(norm, A, U, S, VT );
	bool res = true;
	res &= svd.computaionSuccessed();
	double limit = 1e-12;
	for(int i = 0; i < nNorms; i++){
		res &= norm[i] < limit;
	}
//	cout << "Extra mem " << svd.getMemUsage() << " Bytes" << endl;


	// check the rusult up to some relative accuracy
/* 	double ref_s[] =   {  	 2.7468732418221837e+01, 2.2643185009774690e+01,
	 8.5583882284825830e+00, 5.9857232015121360e+00, 2.0148996587157564e+00 };
	bool res = true;
	const double reltol = 1e-14;
	for (int i = 0; i < n; i++) 
		res &= fabs(Sigma.data[i] - ref_s[i]) < reltol * ref_s[i];
	if (!res) {
	  cout << endl << svd.getName() << " does not passed the test1" << endl;
	}
	return res;
*/
	return res;
}
FullSVDBenchmark::FullSVDBenchmark(FullSVD & svd_):Benchmark(svd_) {
	svd = &svd_;
	set(baseName, "Benchmark FullSVD");
	const char * s1 = baseName;
	const char * s2 = svd->getName();
	char s3[strlen(s1) + strlen(s2) + 2];
	sprintf(s3, "%s %s", s1, s2);
	setName 				( s3);
	setDescription  	("Measures performance of registered svd implementations");
	setTestAuthor		("Yordan Madzhunkov");
	setDeveloperAuthor("Yordan Madzhunkov");
	nRuns = 10;
	nConfigurations = 10;
	setConfigurationEachRun = true;
	generator = new MatrixGenerator();
	generator->setZeroProb  (0.0);
	A = A2 = U = V = S = NULL;
	setCompileTime(__DATE__, __TIME__);	
	setCompileTime(svd->getCompileTime());; 
// METRICS 
	addTimerMetrics();
	addMemoryMetrics();
	addSuccessMetrics();
	for (int i = 0; i < nSvdNorms; i++) svdNormMetrics[i] = 0;
	addSVDNormsMetrics();
}
FullSVDBenchmark::~FullSVDBenchmark() {
	delObj (A );
	delObj (A2);
   delObj (U );
   delObj (S );
   delObj (V );
	delObj (generator);
}
int FullSVDBenchmark::getProblemSize(const int configuration) const {
	const int problemSize[] = { 10, 25, 50, 75, 100, 150, 200, 300, 400, 500, 600, 700, 1000 };
	return problemSize[configuration];

};
void FullSVDBenchmark::setRunConfiguration(const int n) {
	int N = getProblemSize(n);
	if (A ) delete A ; A  = 0;	A  = new Matrix(N, N, General);
	if (A2) delete A2; A2 = 0;	A2 = new Matrix(N, N, General);
	if (U ) delete U ; U  = 0;	U  = new Matrix(N, N, Orthogonal);
	if (U ) U->makeColMajor();
	if (S ) delete S ; S  = 0;	S  = new Matrix(N, N, Diagonal  );
	if (V ) delete V ; V  = 0;	V  = new Matrix(N, N, Orthogonal);
	if (A && U && S && V) {
		generator->generate(*A);
		A2->copyOnlyData(*A);
	}

}
void FullSVDBenchmark::prepareRun() {}
void FullSVDBenchmark::makeRun() {
	if (A && U && S && V && svd) { 
		svd->saveCompute(*A, *U, *S, *V);
	}
}

void FullSVDBenchmark::addSVDNormsMetrics() {
	for (int i = 0; i < nSvdNorms; i++) 
	if (svdNormMetrics[i] == 0) {
		svdNormMetrics[i] = new Metric<double> (nRuns);	
		if (svdNormMetrics[i]) 
			addMetric(*(svdNormMetrics[i]), normNames[i]);
	}
} 

void FullSVDBenchmark::measureSVDNorms() {
	if (A && U && S && V ) {
		for (int i = 0; i < nSvdNorms; i++)
			if(!svdNormMetrics[i]) return;
		double svdNorms[nSvdNorms];
		getResidual(svdNorms[0], *A2, *U, *S, *V);
		double tmp;
		getOrth    (svdNorms[1], *U); 	getOrth    (tmp, *V);
		svdNorms[1] = fmax(tmp, svdNorms[1]);
		for (int i = 0; i < nSvdNorms; i++)
			svdNormMetrics[i]->addMeasurement(svdNorms[i]);
	}

}
void FullSVDBenchmark::finalizeRun() {
	measureSVDNorms();
}


FullSVDScaling::FullSVDScaling(FullSVD & svd_): Benchmark(svd_) {
	svd = &svd_;
	set(baseName, "Scaling FullSVD");
	const char * s1 = baseName;
	const char * s2 = svd->getName();
	char s3[strlen(s1) + strlen(s2) + 2];
	sprintf(s3, "%s %s", s1, s2);
	setName 				( s3);
	setDescription  	("Measures scaling of registered parallel svd implementations");
	setTestAuthor		("Yordan Madzhunkov");
	setDeveloperAuthor("Yordan Madzhunkov");
	nRuns = 10;
	nConfigurations = ThreadPool::getOnlineCores();
	setConfigurationEachRun = true;
	generator = new MatrixGenerator();
	generator->setZeroProb  (0.0);
	A = A2 = U = V = S = NULL;
	int N = 500;
	A  = new Matrix(N, N, General);
	A2 = new Matrix(N, N, General);
	U  = new Matrix(N, N, Orthogonal);
	S  = new Matrix(N, N, Diagonal  );
	V  = new Matrix(N, N, Orthogonal);
	if (U ) U->makeColMajor();
	if (A ) generator->generate(*A);
	if (A2 && A) A2->copyWithData(*A);
	setCompileTime(__DATE__, __TIME__);	
	setCompileTime(svd->getCompileTime());; 
// METRICS 
	addTimerMetrics();
	addMemoryMetrics();
	addSuccessMetrics();
	for (int i = 0; i < nSvdNorms; i++) 
		svdNormMetrics[i] = 0;
	addSVDNormsMetrics();
	xlabel = "Number of cores";
}
FullSVDScaling::~FullSVDScaling() {
	delObj (A );
	delObj (A2);
   delObj (U );
   delObj (S );
   delObj (V );
	delObj (generator);
}
void FullSVDScaling::setRunConfiguration(const int n) {
	int threads = getProblemSize(n);
	if (svd && svd->getCanBeParallel()) {
		svd->setNumberOfThreads(threads);
	};
}
int  FullSVDScaling::getProblemSize(const int configuration) const {
	return configuration + 1;
}
void FullSVDScaling::prepareRun() {
	if (A2 && A) A->copyWithData(*A2);
}
void FullSVDScaling::makeRun() {
	if (A && U && S && V && svd) { 
		svd->saveCompute(*A, *U, *S, *V);
	}
}
void FullSVDScaling::finalizeRun() {
	measureSVDNorms();
}
void FullSVDScaling::addSVDNormsMetrics() {
	for (int i = 0; i < nSvdNorms; i++) 
	if (svdNormMetrics[i] == 0) {
		svdNormMetrics[i] = new Metric<double> (nRuns);	
		if (svdNormMetrics[i]) 
			addMetric(*(svdNormMetrics[i]), normNames[i]);
	}
}
void FullSVDScaling::measureSVDNorms() {
	if (A && U && S && V ) {
		for (int i = 0; i < nSvdNorms; i++)
			if(!svdNormMetrics[i]) return;
		double svdNorms[nSvdNorms];
		getResidual(svdNorms[0], *A2, *U, *S, *V);
		double tmp;
		getOrth    (svdNorms[1], *U); 	getOrth    (tmp, *V);
		svdNorms[1] = fmax(tmp, svdNorms[1]);
		for (int i = 0; i < nSvdNorms; i++)
			svdNormMetrics[i]->addMeasurement(svdNorms[i]);
	}

}

NormSVDTest normSVDTest;
